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Abstract 

Reptation,  or  accommodation,  is  manifested  in  ferromagnetic  materials  in  a  variety  of  operating  regimes  and 
hence  must  be  incorporated  in  models  used  for  comprehensive  material  characterization  or  model-based  control 
design.  Because  the  microscopic  mechanisms  which  cause  reptation  are  complex,  we  characterize  the  effect  in  a 
phenomenological  macroscopic  manner  within  the  context  of  a  homogenized  energy  framework  for  ferromagnetic 
hysteresis.  Attributes  of  the  model  are  illustrated  through  comparison  with  experimental  data. 
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1.  Introduction 

Reptation  in  ferromagnetic  materials  is  the  phenomenon  in  which  repeated  biased  minor  hysteresis  loops 
with  identical  input  fields  yield  output  inductions/magnetizations  that  shift  with  each  loop.  This  is  also  referred 
to  as  accommodation  and  is  shown  for  a  nickel  rod  in  Figure  1.  Each  shift  is  less  than  the  previous  one  so 
that  the  material  approaches  a  steady-state  condition  if  enough  loops  are  performed.  However,  this  steady- 
state  is  temporary  since  upon  traversing  a  major  loop  and  returning  to  the  previous  input  level,  reptation  is 
again  observed.  For  this  reason,  reptation  cannot  be  “trained”  out  of  a  material,  and  must  be  included  in  any 
mathematical  model  when  high  precision  is  needed. 

The  physical  causes  of  reptation  are  not  well  understood.  Although  the  behavior  bears  a  resemblance  to 
thermally-induced  creep  or  after-effects  [2,  3,  4],  it  cannot  be  fully  explained  by  thermal  effects  alone.  It  is  believed 
to  be  related  to  moment  interaction  fields,  but  these  themselves  are  typically  only  described  in  a  phenomenological 
manner.  Hence  reptation  has  been  characterized  in  a  phenomenological  manner  as  a  combination  of  history- 
dependent  and  history-independent  Preisach  models  [8].  Here  we  take  an  alternative  approach  in  which  we 
relate  reptation  to  the  underlying  energy  of  the  material.  This  is  done  within  the  context  of  the  homogenized 
energy  model,  where  energy  principles  on  a  lattice  level  and  statistical  distributions  are  combined  to  provide  a 
macroscopic  model. 

The  homogenized  energy  framework  for  ferromagnetic  hysteresis  is  summarized  in  Section  2  and  extended 
in  Section  3  to  incorporate  reptation.  Attributes  of  the  extended  model  are  illustrated  in  Section  4  where  the 
model  response  is  compared  with  experimental  data  from  a  nickel  rod. 


2.  Homogenized  Energy  Framework 

As  detailed  in  [2,  5,  6,  7],  the  homogenized  energy  model  for  ferromagnetic  materials  is 
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[M(H))(t)=  /  [MiH  +  H^HMv^HMH^dHjdHc  (1) 
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Figure  1:  Magnetic  field  versus  magnetization  for  a  steel  rod.  Both  minor  loops  were  produced  with  identical 
input  field  levels  and  the  shift  in  the  resulting  magnetization  is  due  to  reptation  effects. 


where  H  is  the  input  magnetic  field,  M  is  the  macroscopic  magnetization,  Hc  and  Hi  are  local  coercive  and 
interaction  fields,  and  vc,  vi  are  associated  densities  which  incorporate  the  effects  of  material  nonhomogeneities 
and  variable  effective  fields  ( He  =  H  +  Hi).  The  integrals  are  typically  approximated  using  Gaussian  or  Newton- 
Cotes  quadrature  techniques.  To  illustrate,  let  Hc[i],  i  =  1, . . .  ,NC  and  Hi\j],j  =  1, . . . , iVj  denote  the  abscissas 
(quadrature  points)  and  wc[i],  Wi\j)  denote  the  respective  quadrature  weights.  Approximation  of  the  integrals 
in  (1)  then  yields  the  discretized  model 

Nc  Nj 

[M(H)](t)  =  +  HI[j];Hc[i])](t)i'c(Hc[i])uI(HI[j])wc[i]wI\j] 

*= i  j=i 

(2) 

Na  N: 

=  ^2^2[m(h  +  Hi\j];Hc[i])](t)wc[i]wI\j], 

*= 1  j= 1 

where  wjc[*]  =  Wc[i\vc(Hc[i})  and  wi\j]  =  wi[j]vi{Hi[j\). 

The  kernel  M  is  formulated  by  balancing  the  Gibbs  and  relative  thermal  energies.  For  simplicity,  we  assume 
that  the  operating  temperature  of  the  device  is  approximately  constant.  For  the  one-dimensional  model,  the 
Gibbs  energy  is  then  given  by 

G  =  ip(M)  -  no HM  (3) 

where  Ho  is  the  permeability  of  free  space,  and  the  Helmholtz  energy  is 

f  \r]{M  +  Mr)2,  M  <  -Mr 

=  J  h(M  -  mR )2>  M  >  M/  (4) 

{  i t,(Mi  -  \M\  <  Mj. 

Here  Mu  denotes  the  remanence  magnetization,  M/  is  the  inflection  point  and  r/  is  the  approximate  slope 
after  switching  occurs.  The  Gibbs  energy  is  balanced  with  the  thermal  energy  via  the  Boltzmann  relation 

H(G)  =  Ce~G^M^kT  (5) 

where  V  is  a  reference  volume,  T  is  the  temperature  in  degrees  Kelvin,  k  is  Boltzmann’s  constant,  and  C  is 
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chosen  to  ensure  integration  to  unity.  As  detailed  in  [5],  this  yields  the  average  magnetization  relations 


<M+> 


Me~G(H,M)V/kTdM 


I M0 


/  e~G^H^v/kTdM 
I M0 


,Mn 


Me-G{H,M)V/kTdM 


<M_>  = 


rMn 


-G(H,M)V/kT  dM 


(6) 


for  positively  and  negatively  oriented  moments,  respectively.  In  practice,  the  positive  or  negative  inflection  point 
(±M/)  is  substituted  for  Mq  for  computational  efficiency. 

Let  x+  denote  the  fraction  of  magnetic  moments  oriented  in  the  positive  direction.  Then  X-  =  1  —  x+ 
moments  are  oriented  negatively,  and  the  kernel  or  hysteron  is  given  by 

M  =  x+  (M+)  +  (1  -  x+)  (M_) .  (7) 

The  moment  fraction  x+  evolves  according  to  the  differential  equation 

X+  =  -P+-X+  +P-+{1  ~  X+)  (8) 


where  p. ^ denotes  the  likelihood  of  a  moment  switching  from  positive  to  negative,  and  p j_  is  the  likelihood  of 

switching  from  positive  to  negative.  These  likelihoods  are  quantified  by  the  relations 
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e—G(H,M)V/kT  dM 
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where  T  denotes  the  material-dependent  relaxation  time  (i.e.,  ui  =  ^  quantifies  the  frequency  at  which  moments 
attempt  to  switch)  and  e  is  a  small  positive  constant.  The  value  of  e  is  arbitrary  so  long  as  it  is  significantly 
smaller  than  Pj.  However,  the  values  of  T  and  V  implicitly  depend  on  e.  Thus,  the  integral  does  not  limit  to  0  as 
e  decreases.  As  detailed  in  [5],  point  evaluation  ( e-G(.HMi)v lkT )  jn  the  numerator  is  sometimes  employed  rather 
than  the  integrals  containing  e,  especially  in  higher  dimensions.  This  yields  nearly  identical  model  predictions 
but  must  be  interpreted  as  approximate  from  a  statistical  sense  since  the  likelihood  of  any  point  occurring  within 
a  continuum  is  0. 


3.  Reptation  Model 

3.1.  Theoretical  Model 

The  physical  causes  of  reptation  are  not  well  understood  but  are  known  to  be  dependent  on  the  electron  spins 
of  the  magnetic  moments.  Repeated  minor  loops  over  the  same  input  field  lessen  the  spin,  which  in  turn  changes 
the  level  of  interaction  between  a  given  moment  and  its  neighbors.  Thus,  reptation  is  tied  to  the  behavior  of  the 
interaction  density  zq.  However,  this  does  not  in  itself  give  sufficient  detail  to  implement  a  model.  As  such,  we 
impose  the  empirical  constraints: 

•  reptation  behavior  approaches  a  steady-state  response, 

•  reptation  behavior  resets  so  that  after  a  major  loop,  any  previous  training  is  lost,  and 


•  reptation  behavior  is  dependent  on  input  field  levels;  training  at  one  input  level  has  no  effect  on  reptation 
amounts  at  other  input  levels. 

To  this  end,  reptation  is  considered  by  changing  the  density  vj  in  (1)  and  (2).  Let  vnom  represent  the  nominal 
or  steady-state  interactive  field  density.  This  is  identified  as  a  model  parameter  in  the  same  manner  iq  was 
considered  previously.  The  field-dependent  density  is  then  defined  as 


MH) 


( _ H _ \ 

V(1  +  «)A  J 

(l  +  a)*k 


(10) 
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where  a  >  0  and  A  €  [0, 1)  are  model  parameters  and  k  £  N  is  an  internal  counter.  That  is,  we  define  the 
interaction  field  density  as  a  scaled  multiple  of  the  nominal  interaction  field  density.  This  scaled  version  is 
depicted  for  various  k  values  in  Figure  2  and  in  appearance  is  similar  to  stretching  the  density.  The  denominator 
in  (10)  assures  that  the  area  under  the  density  remains  constant  as  k  changes  —  the  overall  amount  of  influence 
exerted  by  a  moment  on  others  is  constant  but  this  influence  becomes  more  localized  as  reptation  increases.  The 
additional  two  parameters  a  and  A  quantify  the  difference  between  the  nominal  and  beginning  (A  =  1)  densities 
and  the  rate  at  which  isj  approaches  vn0m-  Note  that  other  methods  of  convergence  to  a  nominal  density  are 
possible  —  e.g.,  a  normal  density  whose  variance  decreases  to  some  constant  value  -  but  (10)  is  chosen  here  for 
its  application  to  general  densities  and  its  computational  simplicity. 

The  requirement  that  reptation  be  related  to  moment  interactions  and  that  behavior  approach  a  steady  state 
are  satisfied  in  (10).  The  remaining  two  requirements  effect  the  choice  of  k.  Reptation  behavior  is  dependent 
on  input  field  levels  -  the  behavior  at  one  field  level  is  independent  of  that  given  at  another.  Thus,  k  must  be  a 
function  of  the  input  field  H.  We  let  A:  be  a  counter  over  H: 

k(H)  =  number  of  times  input  level  has  passed  through  a  value  of  H.  (11) 

Note  that  k(H)  >  1  whenever  it  is  employed.  The  evolution  of  (11)  is  depicted  for  an  example  input  in  Figure  3. 

It  is  apparent  that  k(H )  is  discontinuous  which  introduces  a  discontinuity  in  the  density  (10).  However,  it 
important  to  note  that  this  discretization  does  not  introduce  discontinuities  in  the  magnetization.  Consider  for 
example  the  increasing  H  portion  of  a  major  loop  together  with  several  minor  loops.  As  vi  approaches  vnom  the 
effective  field  approaches  the  applied  field.  A  negatively  oriented  moment  may  in  the  process  cross  its  coercive 
point  and  switch  to  positive.  However,  as  long  as  the  applied  field  is  greater  than  the  negative  coercive  point, 
no  positively  oriented  moment  will  switch  to  negative.  Thus,  the  magnetization  shifts  upward  as  k  increases, 
i.e.,  as  minor  loops  are  traversed.  The  only  change  to  the  magnetization  results  from  some  moments  switching 
earlier  than  they  would  have  otherwise.  A  more  interesting  case  occurs  when  leaving  minor  loops.  Here  k  drops 
to  a  lower  value  (presumably  1).  The  interaction  field  increases  when  this  occurs.  No  moments  will  switch  from 
negative  to  positive  at  this  point  —  any  whose  effective  field  is  moved  above  the  coercive  point  in  this  process 
would  have  already  switched  when  entering  the  first  minor  loop  —  but  a  few  moments  that  switched  positive 
may  have  an  effective  field  that  goes  sufficiently  negative  to  cross  the  negative  coercive  point  and  switch  back  to 
negative  orientation.  This  would  result  in  a  magnetization  that  exits  the  minor  loops  at  a  value  below  what  it 
obtained  in  the  minor  loops.  This  prediction  is  observed  in  data  as  can  be  seen  in  Figure  1. 

Behavior  may  be  reset  by  setting  k(H)  =  0  for  all  H  before  the  increment  is  performed.  The  empirical 
constraints  do  not  specify  when  such  a  reset  occurs,  only  that  it  must  occur  at  least  once  per  major  loop. 
However,  the  reset  is  by  nature  an  irreversible  behavior.  The  homogenized  energy  framework  includes  one  other 


H 

Figure  2:  Distribution  vi  given  in  (10)  as  a  function  of  k. 
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Figure  3:  Depiction  of  the  evolution  of  the  counter  k  as  a  function  of  the  field  H. 


irreversible  behavior  —  the  switch  of  a  moment’s  orientation  from  negative  to  positive  or  positive  to  negative. 
We  therefore  group  these  irreversible  behaviors  and  reset  (11)  whenever  the  moment  switches.  This  implies 
Nc  x  Ni  separate  instances  of  (11)  must  be  employed,  one  for  each  combination  of  Hc  and  Hi  quadrature  points. 
In  Section  3.2,  we  introduce  an  approximation  to  lessen  this  requirement. 

3.2.  Implementation  Issues 

The  density  in  (10)  is  typically  not  known  a  priori.  It  may  then  appear  that  quadrature  points  and  weights 
must  be  calculated  for  vi  at  each  timestep.  However,  this  is  not  the  case.  If  Hinom  and  winorn  represent  the 
quadrature  points  and  products  of  quadrature  weights  and  density  heights  (wi  =  Wi  ■  vj)  for  the  nominal  density 
Vnom,  then  the  quadrature  points  and  weights  for  vi  can  be  calculated  as 

Hi  =  Hinorri  (1  +  a)^  ,  Wi  =  u>i„ofn-  (1^) 

Thus,  the  weights  wi  remain  unchanged  from  the  nominal  density  —  the  division  in  (10)  cancels  with  the  change 
in  stepsize  between  quadrature  points.  The  quadrature  points  are  multiples  of  those  used  for  the  nominal  density. 

It  is  impossible  to  store  the  continuum  required  for  (11).  However,  H  is  bounded  by  material  properties, 
allowing  (11)  to  be  discretized  into  an  array  of  integer  values.  For  a  given  input  H,  the  closest  value  in  the 
table  is  used.  The  size  of  the  table  therefore  depends  on  the  desired  accuracy,  but  a  few  hundred  values  to  a 
thousand  values  should  be  sufficient  for  most  applications.  Computational  savings  are  achieved  by  realizing  that 
k  is  typically  small.  If  we  enforce  k  <  255,  for  example,  one  byte  of  storage  can  be  used  for  each  value.  This 
suggests  table  sizes  are  less  than  one  kilobyte,  which  is  easily  attainable  on  most  platforms.  Furthermore,  the 
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multipliers  (1  +  a)A  may  be  computed  a  priori  and  stored  in  a  vector,  leaving  only  an  increment,  a  multiply, 
possibly  a  memory  reset  to  be  performed  for  each  table  at  each  timestep. 

However,  the  theory  still  requires  one  table  of  k  values  for  each  element  in  the  coercive  x  interactive  mesh. 
This  can  be  reduced  by  approximating  this  with  one  table  per  interactive  density  quadrature  point.  This  table 
is  then  shared  across  all  the  coercive  field  quadrature  points  and  reset  whenever  the  largest  quadrature  point  is 
passed  (Recall  that  the  coercive  distribution  is  bounded  by  a  decaying  exponential,  so  the  largest  quadrature 
point  is  finite).  This  will  typically  require  only  15-80  tables  of  k  values,  instead  of  potentially  thousand. 


4.  Model  Attributes  and  Validation 

To  illustrate  attributes  of  the  model,  we  summarize  both  simulation  results  and  a  model  comparison  with 
the  nickel  data  plotted  in  Figure  1.  To  construct  the  model,  parameters  were  initiall  identified  through  a  least 
squares  fit  to  the  data. 

To  illustrate  the  qualitative  behavior  of  the  model,  the  field  shown  in  Figure  4(a)  was  input  to  the  model  to 
generate  the  magnetization  plotted  in  Figure  4(b)-(d).  The  time-dependent  plot  in  Figure  4(c)  and  expanded 


(a)  (b) 


(c)  (d) 

Figure  4:  (a)  Generated  input,  (b)  and  (d)  magnetic  field  versus  modeled  magnetization,  and  (c)  changes  in  the 
modeled  magnetization  as  a  function  of  time. 
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perspective  in  Figure  4(d)  illustrate  the  repration  behavior  that  occurs  during  multiple  minor  loops  traversals 
as  well  as  the  approach  to  a  steady  state  minor  loop.  This  qualitatively  matches  measured  magnetic  behavior 
(e.g.,  see  [1], 

The  measured  field  shown  in  Figure  5(a)  was  subsequently  input  to  the  model  to  yield  the  major  and  minor 
loops  shown  in  Figure  5(d)  and  (e).  Comparison  with  the  experimental  values  in  Figure  5(b)  and  (c)  —  or, 
equivalently,  Figure  1  —  illustrate  that  the  model  is  quantifying  the  measured  magnetization  behavior  of  the 
nickel  sample. 


5.  Concluding  Remarks 

In  this  paper,  we  have  extended  the  modeling  framework  of  [5,  6,  7]  for  characterizing  ferromagnetic  hysteresis 
to  incorporate  repration  or  accommodation  effects.  This  is  necessary  both  for  fundamental  material  characteri¬ 
zation  and  for  quantifying  the  behavior  of  magnetic  transducers  employed  as  actuators  or  sensors  in  high  drive 
regimes.  These  phenomena  were  incorporated  by  modifying  the  interaction  field  density  to  include  the  migra¬ 
tion,  approach  to  steady  state  behavior,  and  resetting  as  minor  loops  are  traversed  and  exited.  Aspects  of  the 
extended  framework  were  illustrated  through  comparison  with  experimental  nickel  data. 
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Figure  5:  (a)  Input  magnetic  field,  (b)  and  (c)  data  measured  from  a  nickel  200  rod,  and  (d)  and  (e)  corresponding 
model  fit. 


